Anomalous orbital magnetism in Dirac-electron systems: 
Role of pseudo-spin paramagnetism 
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The orbital diamagnetic susceptibility is calculated in monolayer and bilayer graphenes with band 
gap as well as in three-dimensional Dirac systems. It is demonstrated that the pseudo-spin degree 
of freedom such as valleys produces paramagnetic susceptibility in an equal manner as the real spin 
dominating over the Landau diamagnetism. The pseudo-spin paramagnetism explains the origin of 
a singular diamagnetism which is present only in the band-gap region and disappears rapidly inside 
the conduction and valence bands. 



I. INTRODUCTION 

The magnetism of conventional metal is composed of 
two different contributions, the spin component known as 
the Pauli paramagnetism and the orbital component as 
the Landau diamagnetism. In condensed matter systems, 
the orbital magnetism sensitively depends on the detail 
of the electronic band structure, and sometimes largely 
deviates from the Landau diamagnetism. Particularly, 
narrow gap materials such as graphite^"— or bismuth^^— 
exhibit a singular behavior in the orbital susceptibility 
near the energy gap. In this paper, we show that the 
anomalous orbital magnetism in narrow gap systems can 
be understood in terms of pseudo-spin paramagnetism, 
which arises from the extra degree of freedom in the or- 
bital motion of electrons. 

Graphene monolayer recently fabricated^Ti^ is a zero- 
gap system in which the conduction and valance bands 
stick together at K and K' points located at inequiv- 
alent corners of the Brillouin zone, called valleys jiiiSrii 
The system is characterized by chiral quasiparticles with 
opposite chirality in each valley and a linear dispersion 
reminiscent of massless Dirac fermions. At the Dirac 
point where two bands cross each other, the magnetic 
susceptibility has a singularity expressed as a delta func- 
tion in Fermi energy ep and disappears otherwise ii^^— 

Bilayer graphene composed of a pair of graphene 
layers^"— has a zero-gap structure with quadratic 
dispersion)^"— leading to a less singular, logarithmic 
peak of the susceptibility] The orbital magnetism 
was also studied for related materials, such as graphite 
intercalation compounds )^i^"— carbon nanotubc)^"— 
few-layer graphenes f^^iiiii^ and organic compounds hav- 
ing Dirac-like spectrum4^ In both monolayer-i2i^ and 
bilayer graphenes;^"— i^^i^i^ certain asymmetric poten- 
tial opens an energy gap at the band touching point. 

In this paper, we calculate the orbital magnetism of 
several Dirac-like systems with gap and show that the 
pseudo-spin degree of freedom such as valleys in graphene 
produces paramagnetism in the same manner as the real 
spin and gives an essential contribution to the singular 
diamagnetic behavior. In Sec. |lTl we calculate the sus- 
ceptibility of the monolayer graphene with varying gap. 



The singular susceptibility-change in varying e^? is un- 
derstood in terms of valley-induced paramagnetism. We 
extend the analysis to the bilayer graphene in Sec. IIIII 
and to a three-dimensional Dirac system corresponding 
to bismutbi"— 1^ in Sec. IIVI A brief conclusion is pre- 
sented in Sec. El 



II. MONOLAYER GRAPHENE 

Graphene is composed of a honeycomb network of car- 
bon atoms, where a unit cell contains a pair of sublattices, 
denoted by A and B. Electronic states in the vicinity of 
K and K' points in the Brillouin zone are well described 
by the effective mass approximation»iii2r2i Let \A) and 
\B) be the Bloch functions at the K point, correspond- 
ing to the A and B sublattices, respectively. In a basis 
{\A), \B)), the Hamiltonian for the monolayer graphene 
around the K point becomesiiiSril 



H^=(_r (1) 



where v is the velocity, 'k± = tTj, ± iny, and tt — 
—ih^ + (e/c)A with vector potential A giving external 
magnetic field B = V x A. In the following, we shall 
completely neglect the spin Zeeman energy because the 
spin splitting is much smaller than Landau-level separa- 
tions. The Hamiltonian at the K' point is obtained by 
exchanging tt± in Eq. ([T]). 

The diagonal terms ±A represent the potential asym- 
metry between A and B sites, which opens an energy gap 
at the Dirac point ji^>^ This can arise when graphene is 
placed on a certain substrate material. In fact, band-gap 
opening is observed in graphene epitaxially grown on a 
SiC substrate.—"^ From a theoretical point of view, the 
singular behavior in ideal graphene with vanishing gap is 
better understood by taking the limit A — >■ 0, as will be 
shown below. We can safely assume A > without loss 
of generality. 

The energy band at B = is given by 
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A2 
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with electron momentum p = {px,Py) andp = ^Jp: 
The density of states ii 
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where gs — 2 and — 2 represent the degrees of freedom 
associated with spin and valley, respectively, and 9{t) is 
a step function, defined by 



m 



1 (i>0); 
{t<0). 



(4) 



The Landau-level spectrum can be found using the re- 
lation 7r+ = {'\/2h/lB)a) and 7r_ = {y/2h/lB)a, where 
Ib = a/ ch/{eB) is magnetic length and and a are rais- 
ing and lowering operators for usual Landau-level wave 
functions, respectively. The eigenenergy at K point be- 
comes 



sgn_(n)V(fiwBFRTA2 (n = 0,±l,±2,---),(5) 



where huJB = \/2hv/lB and 



sgn_t(n) 



+1 (n>0); 
±1 (n = 0); 
-1 {n < 0). 



The corresponding wavefunction is 

,K _ /^sin(Q;„/2)(?!)|„|_] 



where a„ satisfies 



\^ cos(a„/2) 



hujf 



sma^ 



cos ttr 



n\ sgn_(n) 



y/{hLUBr\n\+A^' 
A sgn_ (n) 



(6) 

(7) 

(8) 
(9) 



and 0„ is the usual Landau-level wave function, where 
(/)„ with n < should be regarded as 0. 

The Landau level n = lies just at the top of the 
valence band, i.e., Eq = —A because ao = 0, and its 
amplitude is only at the B site. Similarly, for Landau 
levels lying in the vicinity of the valence-band top, i.e., 
n < satisfying hLjjB\/\n\ <C A, we have a„ « 0, showing 
that the amplitude of the wave function is significant only 
at the B site. For those in the conduction band n > 0, 
on the other hand, we have a„ « tt, showing that the 
amplitude is significant only at the A site. 

For the K' point, on the other hand, the eigenenergy is 
given by Eq. ([5]) with sgn_ (n) being replaced by sgn_|_ (n) 
and the eigenfunction is given by 



K' ^ ( sin(a;,/2) (/)|„| 
l^cosK/2)0|„|_i 



(10) 



where is obtained from Q!„ by replacing sgn_(n) by 
sgn^(n). The Landau level n = lies just at the bottom 



of the conduction band, i.e., = -I- A because a'^ = tt, 
and its amplitude is only at the A site. Similarly, for 
low-lying Landau levels in the conduction band n > 0, 
we have a'^ ~ tt, showing that the amplitude of the wave 
function is significant only at the A site. For those in the 
valence band n < 0, on the other hand, we have a'j ~ 0, 
showing that the amplitude is significant only at the B 
site. 

The Landau levels of n 7^ are doubly degenerate 
between the K and K' valleys, while those of n = are 
not. Therefore, by defining 



A2, 



(11) 
(12) 



the thermodynamical potential at temperature T then 
becomes 



71=0 



where /3 — l/fc^T and 

(^(e)-log[l + e-^(^-^)], 



(13) 



(14) 



with ^ being the chemical potential. 

In weak magnetic field, using the Euler-Maclaurin for- 
mula, the summation in n in Eq. (|13p can be written as 
an integral in continuous variable x and a residual term 

as 



n 



9vgs 



1 



E 



ip[es{x)]dx 



(hujB)'^ dip[es{x)] 



12 dx 
The magnetization is given by 

and the magnetic susceptibility by 



x=0 



OiS'). (15) 



(16) 



lim 

B^o B 



KdB^Jc 



5=0 



(17) 



In Eq. (jl5l) . the first term represents the thermody- 
namic potential in the absence of a magnetic field and 
only the second term depends on the magnetic field. We 
have 



X = I de{ -|^)x(e), 



with 



X(£) = -9vgs 



1 



67rc2 2|A| 



|A|-N), 



(18) 



(19) 



where /(e) is the Fermi distribution function. The sus- 
ceptibility at zero temperature is given by x{^f) with 
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FIG. 1: Orbital susceptibility (solid) and density of states 
(dashed) of monolayer graphene with band gap A. Note that 
the upward direction represents negative (i.e., diamagnetic) 
susceptibility. 



Ep being the Fermi energy. In the limit of A 
susceptibility approaches a delta function 
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0, the 



(20) 



in agreement with the previous result 

The susceptibility and the density of states, given by 
Eqs. ([TS]) and respectively, are shown in Fig.[T] Note 
that the upward direction represents negative (i.e., dia- 
magnetic) susceptibility. The susceptibility is not zero 
at zero electron density, —1 < s/\A\ < +1, because the 
completely filled valence band gives a constant diamag- 
netic susceptibility. When the Fermi energy enters the 
conduction band, the susceptibility jumps downs to zero, 
resulting in zero total magnetism. In the limit of A — >■ 0, 
the step height at the band edge increases in proportion 
to A~^ and the susceptibility approaches a delta func- 
tion. 

Because the Hamiltonian is equivalent to that of a 
Dirac electron with a nonzero mass, the magnetic sus- 
ceptibility around the band edge should correspond to 
that of a conventional electron. This is clearly illustrated 
by the effective Hamiltonian expanded in the vicinity of 
k = 0. For the conduction band, s ~ +1, the effective 
Hamiltonian for the A site near the band bottom (e = A) 
is written apart from the constant energy as 



H 



K 



2K'' 

IK'' 



1 



2m* 



(21) 
(22) 



where /is is the Bohr magneton, given by eh/ (2mc) with 
m being the free electron mass, and we used the relation 
[T^x^T^y] = ifieB/c and defined 



A , 



(23) 



The last term in each Hamiltonian can be regarded as 
the pseudo-spin Zeeman term, where the different valleys 
K and K' serve as pseudo-spin up (^ — +1) and down 
(^ = —1), respectively. This agrees with the Zeeman 
energy expected for an intrinsic magnetic moment, that 
originates from the self-rotation of the wave packet in 
Bloch electroni^i^ The combined Hamiltonian is written 



as 



2m* 2 



(24) 



Obviously, the pseudo-spin Zeeman term gives the 
Pauli paramagnetism and the first term containing tt^ 
gives the Landau diamagnetism in the usual form as 



Xp{e) = {^)%lD{s), 

1 / TTl \ ^ 

Xde) ^ --(-) ADie). 
with density of states 



D{e) 



e{e) 



(25) 
(26) 

(27) 



The total susceptibility xp + XL actually agrees with the 
amount of the jump at the conduction band bottom in x 
of Eq. (fT9| . Because g = 2m/m* in the present case, we 
have XL = — Xp/3 oc 1/m* as in the free electron, giving 
the paramagnetic susceptibility in total. Therefore the 
susceptibility exhibits a discrete jump toward the para- 
magnetic direction when the Fermi energy moves off the 
Dirac point. The jump height goes to infinity as the gap 
closes, because the susceptibility is inversely proportional 
to the effective mass. 

In the original Hamiltonian, the Landau-level energies 
in Eq. ([5]) can be rewritten as 



(n' =0,1,2, 



A2 
••)• 



(28) 



Figure [5] shows energy levels for Tlujb = 2A and the re- 
lationship between the different labeling schemes of Eqs. 
([5]) and (f28|l . For the conduction band, the levels of the 
same n' with opposite pseudo-spins ^ = ±1 share the 
same Landau level function 0„' on the A site, on which 
the states near the conduction-band bottom {e = A) have 
most of the amplitude as has been discussed above. For 
the valence band, similarly, n' describes the index of the 
Landau-level function at the B site. 
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FIG. 2: (color online) Landau-level energies of gapped mono- 
layer graphene with for hujB = 2A. Dashed lines connecting 
the levels of K and K' represent corresponding levels with 
opposite pseudo-spins. 



III. BILAYER GRAPHENE 

Bilayer graphene is a pair of graphene layers arranged 
in AB (Bernal) stacking and includes Ai and Bi atoms on 
layer 1 and A2 and B2 on layer 2^'^— The low energy 
states are again given by the states around K and K' 
points in the Brillouin zone. The Hamiltonian at the K 
point for the basis {\Ai), \Bi), IA2), I-B2)) is given by 



A 







WTT-)- A 71 
71 —A VTT- 

W7r+ -A/ 



(29) 



where A describes potential difference between layer 1 
and 2 (not A and B sites) and 71 represents interlayer 
coupling between Bi and A2 i^^i^Oi^i Hamiltonian at 
the K' point is obtained by exchanging 7r± in Eq. (P^ . 
The energy band at B = is given by 



^+i;V(7? + 4A2) 



l/2\ 1/2 



(30) 



with ^ = zbl.'^'^ The index jj. — +1 and —1 give a pair 
of bands further and closer to zero energy, respectively, 
and s = -1-1 and —1 in each pair represent the electron 
and hole branches, respectively. The band-edge energies 
corresponding to p = are given by jej = e± for /j, = ±1, 
where 



/72 + A2, e- = |A| 
For fi = —1, the band minimum becomes 

7i|A| 



^0 = 



VtF+Ia^ 



(31) 



(32) 
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FIG. 3: (color online) (a) Band structure, (b) density of 
states, and (c) susceptibility of bilayer graphenes with the 
asymmetry gap A/71 = 0, 0.2, and 0.5. Vertical lines in (a) 
indicate the energies of eo, £- and e+ for A/71 — 0.5. The 
upward direction represents negative (i.e., diamagnetic) sus- 
ceptibility in (c). 



which corresponds to an off-center momentum. '^'^ The 
density of states diverges here as D{e) (x {e — Sq)"^^^. 
The energy bands and the density of states with several 
A's are plotted in Figs. |3] (a) and (b), respectively. Ver- 
tical lines in (a) indicate the energies of Eq, and e+ 
for A = 0.571. 

In a magnetic field, the eigenfunction of the Hamil- 
tonian at the K point is written as (ci0„_i, C20n, C30„, 
C4(f)n+i) with integer n > —1. For n > 1, the Hamilto- 



5 



nian matrix for (01,02,03,04) become o^"'^^ 



fc, . /rr A 



K(^ = +) 



71 
-A 



-A 



/ 

(33) 

For n — 0, the first component does not actually exist 
because 0_i = 0. The matrix for (02,03,04) becomes 



/ A hjjB\/n 
hijjB\pn A 

71 
hwB^n + 1 



— 

tif) - 



A 71 
71 —A hjjjB 
hbJB -A 



(34) 



o 




J (+,+) 



1 = 3 

(+■-) 



1 = 2 

(-,-) 



1 = 1 

(-,+) 



For n = — 1, only the component 04 survives and the 
Hamiltonian is 



ii\ = -A. 



(35) 



For the K' point, the eigenfunction is written as 
(ci(?!)„+i,C2(/)„,O30„,C4(/)„_i). For n> 1, the Hamilto- 
nian matrix for (01,02,03,04) is 



J?. 



K' 



( A 

hijJB\fn 



V 



hujBvn 
T A 

71 




1 



For n = 0, the matrix for (01,02,03) becomes 



> 

71 
—A hcuBy/n 
hujB\/n —A J 
(36) 



H, 



K' 



A hujB 
hu!B A 71 

71 -Ay 



and for n = —1, that for ci is 



(37) 



(38) 



If we extend the definition of the matrix of Eq. ([55)1 
to n = 0, its three eigenvalues agree with those of 
and the rest with that of . Similarly, the matrix 
of Eq. (p6)) with n = gives eigenvalues of and 
H^^. Thus we can use Eqs. ^ and §6^ with n > to 
produce the full spectrum. By introducing the pseudo- 
spin variable ^ = ±1, the Hamiltonian is combined into 
a single expression. 







A y^Xn + 

-A 



(39) 



with 



and 



1,. 



(40) 



(41) 



FIG. 4: (color online) Landau-level energies of bilayer 
graphene given by Eq. H39p with A = O.271 and fkos = O.571. 
The dashed horizontal lines represent the Landau level which 
originally belongs to n = — 1 at opposite valley. The quantum 
number (s,/i) is indicated below j. 



We write the eigenvalues of Hf^ as 

ej{xn,^S) (j = 1,2,3,4), 



(42) 



in the ascending order in energy (j = 1 and 2 for valence 
bands and j = 3 and 4 for the conduction bands). The 
second argument in Sj{xn,^S) represents the dependence 
on B which are not included in 

Figure S] shows the example of the Landau- level spec- 
trum at A/71 = 0.2 and Hlvb/ji = 0.5, where the thick 
dashed lines represent the Landau level which originally 
belongs to n = — 1 at opposite valleys. The correspon- 
dence between quantum numbers j and (s,/i) are indi- 
cated in the figure. 

The thermodynamic potential becomes 



n 



9t 



1 



I3 2ttI 



2 X^H^t'^J'^^"'^^)] 



^ i,3 L 

52 d^\e,{xM 



93 [ej (a;, ^5)] da; 



24 



dx 



(43) 



where we used the Euler-Maclaurin formula in the sec- 
ond equation. The first term in the bracket can be trans- 
formed by changing the integral variable from x to e as 

- ^[e,{x,^6)]dx= f{e)nj{e,^6)de, (44) 
P Ja J-00 

where we used ip'{e) = —l3f{e) and defined 



nj{e,^6) = Sj{e,^S)xj{e,^6), 



(45) 



where Xj{e,(,S) is a real and positive solution of e 

£j(x, ^S) and 



Sjie.^S) = sgn(^ 



dxj{e, ^5) 



(46) 
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If there are more than one solution of Xj, we regard nj discrete jump at e+ as 



as their sum. The quantity nj{e,^6)/{4:Trh v ) represents 
the electron density below e for the conduction band and 
the hole density above e for the valence band. 
By expanding 

nj{e,^S) = (e) + nW(e) ^5+^nf\e) + ■ ■ ■ , (47) 

we can further expand of Eq. (|43|) in terms of 6 (x B. 
We have 



X(e) 



2 2 

e w 
9s9v — 5" 



E 



nf\e')de' 



\e — e 



:(0,0)] 



dej{x, 0) 



dx 



x=0 



(48) 



For the Hamiltonian of Eq. (p9| . the eigenequation 
det(e — = can be solved for x (= a;„) as 



x±^e^ + A'± -y (4eA - ^S)"^ + Ajf{e^ - A2), (49) 

which gives Xj{e^^S) when being real and positive. Let 
us first consider the case e > £+, where two conduction 
bands are occupied by electrons. In this case x± are both 
real and positive and we have xi = X2 = 0, a;3 = x-^-, and 
X4 = X-. Then, we have 



2(e2 



A^) 



(50) 



independent of S^S. Therefore, 'J2j (s) identically van- 
ishes, resulting in susceptibility independent of energy in 
the region e > e+. The same is true for e < —e+. Be- 
cause X = for £ = ±oo, i.e., in the case of empty or 
filled band, we can conclude that the susceptibility van- 
ishes for e > and e < — e+ independent of interlayer 
interaction 71 and asymmetry A. 

Similarly, the density of states for |e| > e+ is indepen- 
dent of 71 and A and becomes twice as large as that of 
monolayer. In fact, we have 



d 



^(^)^ a^E"^ (^'0) = 4^■ 



(5l) 



This feature of the density of states is apparent in Fig. [3] 
(b). 

In the vicinity of the bottom of the excited conduction 
band, £ = £+, we have 



,(2) 



is) = 



d^x. 



(5=0 



where we used X-6{x-) = and X-6'{x-) = —S{x_). 
Using Eq. (HSl) . we find that the susceptibility makes a 



X(e4 



gvgs 



0)-x{e 
.2„,2 .AVA2 



0) 



-C2 V 



7? 



2A2 + 72 
7?(2A2+72) 1272 yA^ 



-7f 



(53) 



where the first term in the bracket comes from the inte- 
gral of the delta function in Eq. (|52]) and the second term 
from the step function in Eq. 

Near e+, the eigenstates are given primarily by the 
dimer states composed of \Bi) and IA2). The effective 
Hamiltonian is described by the second order in inter- 
band interation with the conduction-band bottom \Ai) 
and the valence-band top \B2), where each process gives 
a term c>c 7r_|_7r_ or cx 7r_7r-|_. In symmetric bilayer with 
A = 0, the terms tt+tt- and tt-tt-^ have the same coeffi- 
cient and the pseudo-spin Zeeman term identically van- 
ishes. When A becomes nonzero, the two coefficients 
shift from each other linearly in A because of the band- 
gap opening, leading to a nonzero Zeeman term. The 
resulting effective Hamiltonian is given by Eq. (f24| with 



m 



7^2 v/A2 4-72 ^ ^ 4Av/A2-t-72 m 
2w2(2A2h-72)' ^ 2A2-H72 m* 



(54) 



The susceptibility is written as Pauli and Landau mag- 
netism in Eqs. (^5]) and I^E^, respectively, which together 
give a susceptibility jump of Eq. ()53p . The paramagnetic 
component xp is zero at A = and monotonically in- 
creases as A becomes larger. At g* — {2/\/3){m/m*) or 
A « 0.3471, Xp exceeds xl and the susceptibfiity step 
changes from diamagnetic to paramagnetic. In the limit 
A — > 00, we have g* — 2m/ m* as in the monolayer. This 
is to be expected, because the bilayer graphene in this 
limit can be regarded as a pair of independent monolayer 
graphenes, where interlayer coupling 71 opens an energy 
gap at each Dirac point. Similar argument also applies 
to the behavior around 

In the energy region — e_ < e < —Sq near the top of 
the valence band, both xj^. and X- are real and positive, 
giving the states at outer and inner equi-energy circle of 
the band j = 2, respectively. Then we have 



.(2) 



ie) 



— i-x^+x^) 
7?(A2-, 



(5=0 
2^ 



2[(4A2+72)(£2_,2)p/ 



(55) 



When the energy approaches to — £0 from the negative 
(2) 

side, the integral of (e), thus the susceptibility, di- 
verges in positive direction as oc {e + So)~^^'^ in the same 
manner as the density of states. The same divergence oc- 
curs at the bottom of the conduction band, -l-eo, because 
of the electron-hole symmetry. 

Full analytic expression of the susceptibility xi^) is 
complicated and presented in Appendix [Xj Figure [3] (c) 
plots the susceptibility for A = 0, 0.2, and 0.5. In ac- 
cordance with the above analytic consideration, we actu- 
ally observe that the susceptibility vanishes in the region 



7 



£ > e+ and e < — e+ and that the susceptibility step 
at e = e_|_ changes from diamagnetic to paramagnetic 
with increasing A. We also see that the susceptibility for 
A 7^ diverges in the paramagnetic direction at e = ±eo- 



IV. THREE DIMENSIONAL DIRAC SYSTEM 

The results in monolayer graphene in Sec. |lT] can be 
directly extended to three-dimensional Dirac Hamilto- 
nian, which is also known to describe the approximate 
electronic structure of bismuth with strong spin-orbit 
interactioui^"— In bisumuth, electronic states near the 
Fermi level is approximately described by a (4, 4) matrix 
Hamiltonian, given by 



1.5 



n = 







VTTz VTT^ 

A W7r_|- —VKz 

VTTz VT^- —A 

vTr+ —VTT, —A 



(56) 



where four components consist of two orbital and two 
spin degerees of freedom. The density of states at zero 
magnetic field is 



D{s) 



9v9s 



A2 



(57) 



where is the valley degeneracy allowing the presence 
of different k points described by the above Hamiltonian 
in the first Brillouin zone. 

The Landau levels in a uniform magnetic field in z 
direction are given by 
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FIG. 5: Orbital susceptibility and density of states of three- 
dimensional Dirac electron. 



where D{e) = (gs5„/47r2)(2TO*/ft2)3/2^ ^-^^^i m* = 
A/u^. This is nothing but the magnetic susceptibil- 
ity, dominated by the Pauli paramagnetism, of a three- 
dimensional metal with mass m* and g factor g* — 
2m /m* . Figure[5]shows the susceptibility and the density 
of states in the present system. The singular decrease of 
the susceptibility at the band edges is fully understood 
in terms of the appearance of the dominant spin param- 
agnetism inside the band. 



A2 



(n = 0,l,2,...), (58) 

with hbJB = \/2hv/lB, s ~ ±1, and a — ±1. This is 
equivalent to the two-dimensional Dirac system, Eq. (1281) , 
when the term A^ is replaced with A^ -|- u^p^- The sus- 
ceptibility x(£) is thus calculated by integrating Eq. 

in "Pz as 



X(£) 



2Trh 



v'^p^ 



log 



2£c 



A2 



(kl<|A|); 
(kl>|A|), 



(59) 



where Ec is a cut-off energy. In the limit of A — >■ 0, the 
susceptibility at zero energy logarithmically diverges. 

At an energy e just above the band bottom |A|, we 
obtain the paramagnetic contribution 



2 / m 



(60) 



V. CONCLUSION 

We have shown that the orbital magnetism singularly 
dependent on the Fermi level, appearing in narrow gap 
electronic systems described by the Dirac Hamiltonian, 
can be understood in terms of pseudo-spin Pauli para- 
magnetism induced by extra degree of freedom in the 
orbital motion. This has been demonstrated by explicit 
calculations of orbital susceptibility in monolayer and bi- 
layer graphenes with band gap and also in three dimen- 
sional Dirac systems such as bismuth. 
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Appendix A: Susceptibility of bilayer graphene 



and 



Using Eqs. (gH]) and (gH), the susceptibility of x(£) of 
bilayer graphene with energy gap is calculated as 



with 

where 

Fie) 



2 2 

e w 



XO + X-+X+ (kl < eo), 

2F{e) + X-+X+ (eo < kl < 

F{e)+x+ (e-<|e|<£+) 

(£+ < 



7lA2|£| 



(7f + 4A2)/^ 



7? 



4 (7^ + 4A2) 3/2 



4A2e2 



log 2|e|(72 + 4A2) 



(Al) 



(A2) 



Xo 

X- 

x+ 



7?log[472A2(72 + 4A2)] 

4(72 +4A2)3/2 

071 



(A4) 
(A5) 



A2 7A^ 2A2+7? 
'^^'"^ " 7i(2A2+7?) ^ 127xv/A^+7r ^ ^ 



In the symmetric bilayer with vanishing gap, A = 0, in 
particular, we simply get 



Xie) = 5.5.^^(71 - \e\)(\ log ^ + -^.) , (A7) 
7rc^7i V4 71 12/ 



2W72 + A/S?J^l (e2 — A2) + 4A2£2 ^ (A3) which agrees with the previous resultsi^i^S 
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